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Abstract 

We consider a random aggregate of identical frictionless elastic 
spheres that has first been subjected to an isotropic compression and 
then sheared. We assume that the average strain provides a good de- 
scription of how stress is built up in the initial isotropic compression. 
However, when calculating the increment in the displacement between 
a typical pair of contaction particles due to the shearing, we employ 
force equilibrium for the particles of the pair, assuming that the aver- 
age strain provides a good approximation for their interactions with 
their neighbors. The incorporation of these additional degrees of free- 
dom in the displacement of a typical pair relaxes the system, leading 
to a decrease in the effective moduli of the aggregate. The introduc- 
tion of simple models for the statistics of the ordinary and conditional 
averages contributes an additional decrease in moduli. The result- 
ing value of the shear modulus is in far better agreement with that 
measured in numerical simulations. 
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1 Introduction 



Digby (1981) and Walton (1987) considered a random aggregate of frictional 
spheres in which the distribution of contacts was isotropic. They consid- 
ered a random aggregate of identical spheres that was first compressed by 
an average pressure p. They assumed that the relative displacement of the 
centers of two contacting particles was given by the average strain, and they 
obtained expressions for the effective shear modulus fi E and Lame coefficient 
X E . Their expressions for these moduli are 
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where k is the average number of contacts per particle (the coordination 
number) and v is the solid volume fraction. The parameter a describes 
the strength of the transverse stiffness of the grain-to-grain contact; a = is 
appropriate to frictionless interactions (perfect slip), whereas a = 1 describes 
the fully frictional interactions (perfect stick). The effective bulk modulus, 
k e , and the effective Poisson ratio, v E , are given in terms of these by 
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respectively. Note that the bulk modulus, k e , does not depend upon a 
because the transverse forces do not enter at all into this average strain 
approximation. 

The effective moduli of the corresponding aggregate of frictionless spheres 
can be obtained simply from equations (1) and (2) by setting a = 0: 
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The equality of the coefficients is consistent with Cauchy's use of the average 
strain assumption to obtain a single independent modulus for random arrays 
of grains that interact through central forces (e.g., Love (1927), Note B). 
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Jenkins, et al. (1989) compared the predicted values of the effective shear 
and bulk moduli of frictional spheres with the results of computer simula- 
tions and physical experiments on a binary mixture of glass spheres with 
rather large differences in diameters that were isotropically compressed to an 
average pressure of 138 kPa. They found that the effective shear and bulk 
modulus predicted by Digby (1981) and Walton (1987) were, respectively, 
three times and one and one-half times greater than the values measured in 
the experiments and the simulations. Makse, et al. (1999) also compared 
values of the effective moduli with the results of computer simulations. For 
a binary mixture of frictional spheres that differed little in diameter, they 
found the effective shear modulus to be about two-thirds of the value pre- 
dicted using the average strain assumption. They explained this difference 
as being due to the relaxation of the particles associated with their achieving 
equilibrium in the numerical simulation. More surprisingly, they found that 
the effective shear modulus for frictionless spheres in the simulation was less 
than ten per cent of the predicted effective medium value, equation (5). It is 
the central goal of the present article to understand why this is so. (By con- 
trast the bulk modulus agreed reasonably well with equation (3) regardless 
whether there was perfect slip or perfect stick.) 

Because the difficulty with the shear modulus was shown to be due to 
the relaxation of the particles from the average strain, we first perform the 
simplest investigation that allows for some relaxation. From the simulations, 
we know the rest positions of each of the particles, as well as the vectors be- 
tween the centers. Consider a specific particle. We make the approximation 
that when a macroscopic strain increment is applied, the particles in contact 
with it move according to the average strain. Because the specific particle 
is not, in general, in a symmetric environment, it experiences an unbalanced 
force. Consequently, it will, in general, move to a position different than 
that expected from the average strain, in order to reduce the net force to 
zero. So, for the specific particle, we calculate its new position. We next 
calculate the energy stored within each of the contact "springs" for each of 
the particles in the simulation and calculate the total stored energy due to 
the applied strain. We set this equal to the usual expression for strain energy 
and deduce new estimates for the bulk and shear moduli of the aggregate. 
This procedure is detailed in Appendix A. 

It is obvious that such a procedure can only reduce the moduli relative 
to those of the average strain prediction. Note that if we were to neglect 
the relaxation of each particle and assume each that particle sees the same 
coordination number, k, distributed uniformly around it, we would reproduce 
the average strain predictions, equations (1) and (3), as detailed in Norris 
and Johnson (1997). We find that for a static confining pressure of 100 kPa, 
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there is a small reduction of the bulk modulus from 223 MPa predicted by 
the average strain analysis to 206 MPa. There is a much larger reduction 
of the shear modulus from 134 MPa to 100 MPa; however, the results of 
the simulations for the shear modulus give the much, much smaller value 
\i E = 8 ± 3 MPa. We see that relaxation effects at the single particle level, 
while significant, are by no means sufficient to explain the effect. 

We are thus led to consider a more sophisticated theory in which we 
explicitly account for fluctuations in pairs of contacting particles. Here, we 
specialize specifically to the frictionless case, where the reduction in shear 
modulus is most dramatic and for which we can derive an analytic result using 
some fairly weak assumptions. We employ the relatively simple model of a 
static aggregate introduced by Jenkins (1997) in which the assumption that 
the increment in the contact displacement is given by the increment of the 
average strain is relaxed. Instead, the centers of a typical pair of contacting 
particles are assumed to be able to translate in order to equilibrate force, 
while their surrounding neighbors are constrained to move with the increment 
in the average strain. Incremental strains are employed because the contact 
forces are nonlinear functions of the displacement. When the equilibrium 
equations are phrased in terms of increments in displacements, they provide 
linear equations for their determination in terms of the increment in average 
strain. We obtain an approximate analytical solution to the equilibrium 
equations that determine the increments in displacements of the pair in terms 
of the increment of average strain. This solution contains quantities that 
involve the geometry and interactions of a particle with its neighbors; we 
provide relatively simple statistical models for the averages of these and their 
correlations. This permits the calculation of the increment of contact force 
between the pair. Summation of the increments of contact force over all pair 
orientations provides the relationship between the increment of stress and 
the increment in strain and, hence, the effective moduli. In this way, we 
calculate a decrease of the effective shear modulus of about seventy per cent 
from the value predicted by the more elementary theory. 



2 Theory 

We focus our attention on a pair of contacting spheres, label them A and 
B, and denote the vector from the center of A to the center of B by &^ BA \ 
We write the increment F^^in the contact force exerted by particle B on 
particle A in terms of the increment vS BA ^'m the relative displacement of the 
points of contact: 

p(BA )=Ki BA )il (BA) i (6) 
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where K.( BA ^> is the contact stiffness. 

Here, we assume that the contact stiffness is given in terms of the unit 
vector cr BA ) in the direction of by 

K\f A) =K^^df A \ (7) 

where K^ A ^ is the normal contact stiffness, given in terms of the contact 
displacement u^ BA ^ by 

with 

Using the Hertz contact law, the normal displacement can be related to 
average pressure p through the average strain assumption (Jenkins, et al., 
1989) by 
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where d is the sphere diameter. Because computer simulations (Jenkins, et 
al., 1989) indicate that the bulk modulus is rather well predicted by the 
average strain theory, we believe it appropriate to use this expression to 
relate the modulus K N to the pressure p in the initial isotropic state. 

The increment xv^^m. contact displacement may be written in terms 
of the increments and c^in the translations of the centers of the two 
spheres by 

,(BA) _ .(B) _ .(A) 

Alternatively, the relative displacement of the two contacting points may be 
written in terms of the increments in the averages of quantities and their 
fluctuations as 

uf A ^E i3 df A ^5^-8^ A \ (10) 

where E is the increment in the average strain of the aggregate and, for exam- 
ple, Sc^ is the departure of the displacement of sphere B from the average 
strain assumption. The relative displacement can be written more compactly 
by introducing ^ BA ^> = 5c^ — 5c^ A \ the increment in the difference of the 
fluctuations in displacement. 

Given F, the increment T in the stress may be written as the average 
over all N particles in a region of relatively homogeneous strain as 

\ n=l I A=l n=l 
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where is the number of particles in contact with particle A and is 
the volume of the Voroni polyhedron associated with particle A. A contin- 
uous form of this may be written in terms of a distribution function /(d), 
defined so that f(d)dQ is the number of contacts in an element dQ of solid 
angle centered at d : 

tl > = \^ J j[/(d)^(d)didn, 

where the factor multiplying the integral is the half the number of particles 
per unit volume, expressed in terms of the solid volume fraction v. For an 
isotropic distribution of contacts, the distribution can be expressed in terms 
of the coordination number: 




Given the increments E in average strain, the calculation of the increment 
in stress requires the increment in contact force, 

piBA) = tfBA^BA^BA) ^ Jm d^ + Af , 

determined in terms of the fluctuations A for all pairs of particles in the 
region. These fluctuations can be obtained obtained as solutions of the equa- 
tions of balance of force for each of the iV particles. 

3 Pair Fluctuations 

Here, we analyze a far simpler situation in which two contacting particles, 
A and B, have sufficient translational freedom to satisfy force equilibrium. 
In order that the equilibrium equations for the two particles determine these 
translations, we assume that the other particles in contact with the pair 
translate with the average deformation. 

We denote the increment in the translation of center of the n th neighbor 
of particle A by c^. Then, as before, 

AnA) _ (n) _ .(A) 
u i — c i c i ■ 

For n 7^ B, only the fluctuations in the translation of particle A occur; so, 
for these pairs, we may write 

■(«) -(A) _ p AnA) r.(A) _ p AnA) 1 ABA) 1 -y(BA) 
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where E( BA ) = + 5c^ is the increment in the sum of the fluctuations 
in displacement. 

The equations of force equilibrium for particle A are, then, 

= K { N BA) d\ BA) df A) (E jk df A) + Af>) (12) 

NiA) f 11^ 

+ E K^d^d^ [E ]k dt A) + \Kf A ^ - \tf A ^ 

The corresponding equilibrium equations for particle B are obtained by in- 
terchanging A and B, keeping in mind that S AB ^ = — d( BA \ 

The equilibrium equations for the particles A and B lead to a system 
of equations that we use to evaluate the unknown incremental fluctuations 
A( BA } and for the pair BA. In order to phrase the equilibrium equations 

in terms of the neighbors of the individual particles of the pair, we write, for 
example, 

N(A) n (A) 

E K^ A) df A) df A) d { r A) = E K^ A) df A) df A) dt A) -K^df A ^df A hf A) 

n^B n=l 

We then characterize the neighborhoods of particles A and B through the 
tensors 
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The definitions of PsS BA \ A^ AB \ j( BA ^ and J^- 5 ) are based on the existence of 
the contact between A and B and it is this that gives them their directionality. 
Equation (12) can be written in terms of these tensors as 

djg A) E jk +\A^ (Af A) - tf A) Y\K { N BA) dT A) df A) (Af A) + if A) ) = 0. 

Upon interchanging A and £?, an equivalent expression is obtained for the 
force equilibrium for particle B : 

E jk -\^ (Af a) + if A) )-^r ) ^ ) ^r ) (Af a) - tf A) ) = o. 



The summations for A^ BA ^ and A^ AB ^> involve a number of contacts equal, 
on average, to the coordination number; so, up to an error of 1/k, the last 
term in these equations may be neglected. In this case, their solution is 
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Given the tensors A.( BA \ j( BA \ A.( AB \ and 3 {AB \ equations (13) and (14) 
provide the increments in the sum and difference of the fluctuations of the 
pair BA in terms of the increment in average strain. 

Equations (13) and (14) apply to each pair of contacting particles in the 
assembly with orientation near d^ BA \ However, detailed information regard- 
ing the tensors A" 1 J for each such pair is available only from the numerical 
simulations. Consequently, we introduce an average of such tensors. The 
average is taken over all pairs of particles with their orientation in an incre- 
ment of solid angle centered on the unit vector that points in the direction 
indicated by the superscripts. For example, with AQ( BA ^ the increment of 
solid angle centered on the unit vector dS BA ^: 
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where M is the number of pairs in the increment of solid angle. Then, the 
existence of the contact between particle A and B provides a symmetry about 
the plane perpendicular to d^ BA ^: 



f a(AB)\ 1 j(AB) = _( A (BA)\ 1 j(BA) 
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With this, the average value of the sum and difference of the fluctuations 
over all pairs with orientation near d^ BA ^ are 
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In what follows, we also employ the corresponding averages for the indi- 
vidual tensors 
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Becau se of the ro le pl ayed b y cr BA ) in the definition of the average, A( AB ) = 
A( BA ) and J( AB ) = —J( BA ); while the average of 3^ over all orientations 
of the pair is zero. 

In order to evaluate the average on the right hand side of (15), we first 
express the tensors involved as the sum of an average and a fluctuation. For 
example, 

A {BA) =~A^+A (BA) ' 
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Then, in (15), we approximate the inverse of A^ BA ^ in terms of the inverse 
of AMand the fluctuation A^' by 
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We next introduce simple assumptions regarding the distribution of contacts 
and calculate these averages. 



4 Averages 

We introduce an orthogonal Cartesian system with its center coincident with 
the center of the sphere A and characterize a typical contact vector a through 
its equatorial and polar angles (f> and 9: a = (sin # cos 0, sin 6* sin 0, cos 6* ). 
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In calculating the tensors A( BA ) and J( BA ), we replace the summation over 
the discrete contacts of a particle by integration over a contact distribution 
function g( BA \a), defined so that g^ BA \a)dVt is the number of contacts in 
an element of solid angle dVt = sin Odd dip centered at a, given that there is a 
contact at d( BA ) along the polar axis. Then, for example, 



A {BA) 



Kj^df^ + K N 



dQ a , 



where the integration is over all solid angle Q a consistent with the presence of 
particle B. The first term is associated with the presence of this contact. We 
assume that the distribution is uniform in the upper and lower hemispheres. 
That is, in the upper hemisphere, there are, on average, k/2 — 1 contacting 
particles uniformly distributed over the orientations not excluded by the solid 
angle of 7r associated with particle B at the pole; while there are, on average, 
k/2 contacting particles uniformly distributed over the lower hemisphere. 
Then 
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It is a straightforward calculation to determine that 

where an = (19k - 22)/48 and a 2 = (22 - 3Jfc)/16. 
In a similar way, we determine that 
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where wi = (166 - llfc)/128 and cj 2 = -(Jfe + 14)/128. 
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For simplicity, we approximate the tensor A( BA "> by its isotropic part: 
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where ip = k/3. Then its inverse is simply 

-l 



MP =W„)-%. (19) 



With this, the terms in equation (16) that remain to be evaluated are A^^'J^)' 
and A^BWJjBAy. 
Now, by definition, 
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In order to calculate A( BA ) J( BA ~> we introduce the joint probability density 
function F(a,j3), defined so that F(a, f3)dQ, a dQ,p is the fraction of contacts 
with a in dVt a and (3 in dQp, given that there is the contact at d^ BA \ including 
the possibility that the two directions can coincide. Then 
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where the integrals are taken over all solid angles consistent with the presence 
of particle B. Then 
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in which the contributions associated with the presence of the contact at 
<T Byl ) have canceled. 

The unconditional averages are 



Jn f3 
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where uj\ = (38 — Ilk) /128 and uj 2 is defined in (18). 

The joint probability density can be expressed as the product of the simple 
probability g(j3) and the conditional joint probability h a \p(a, j3), defined so 
that h a \p(a, (3)dQ a is the fraction of contacts with a in dQ a , given that j3 is 
in dVLp. 

F(a,P)=g(P)h a \ P (a,P). 

Here h a \p{a, f3) includes the possibility that a equals (3. We take this possibil- 
ity into account explicitly and introduce the conditional probability z a \p{a, f3)dVL a 
that expresses the fraction of contacts at a in dQ a , given that j3 is in dVLp 
with a^p. Then 

/ / F(a,(3)(3 j (3 i a i a m a n dVt a dVti 3 
J dp J d a 

= / g(P)Pjp m p n d^f3+ / / g(p)z a \ l3 (a,p)p j p i a i a m a n dn a dn f3 

J dp J dp J d a 

(22) 

The conditional probability z a \p{a,0) is not uniform; it depends on the lo- 
cations of spheres at (3 and dS BA \ 
^From the above results, we have 

a (ba), a (ba), = K 2 n f f F ^p)p.p 8asaid n p dn a 

J dp J d a 
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with 
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The integrals in (22) and (24) that involve the conditional probability are 
complicated because the limits of integration for 9 a and <p a can depend upon 
Op and <f>p. We illustrate this in the calculation of A( Bj4 )'J( SA )'with the order 
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of integration taken to be d(j) a , d(j)p, dO a , d6p. Then 
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The angle <&(9 a , Op) follows from the fact that the arc on the unit sphere that 
links the centers of the spheres at a and (3 has length d(a, (5) = cos" 1 (a • 0). 
When a is taken to be in the x — z plane, <p a = and <pp = $; then, because 
the spheres are in contact, d = n/3 and 



$ = arccos 



2 cos a cos 0r z 



2 sin a sin Op 

The distribution function z a \p(a,/3) is determined over the range of polar 
angles in each of these integrals in Appendix B. Then, in Appendix C, it is 
used to complete the calculation of A( BA )'J( Bj4 )' and A( Bj4 )'A( Bj4 )'. 
The final results, obtained in Appendix C, are 
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where S^ BA ^ and are given as functions of k and d^ BA ^ in equations 

(41) and (42), respectively, of Appendix C. 

Equations (26) and (27) may be written more compactly as 
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and 

A^A™ = '<% + vJr^r) ■ (29) 

where, with the information provided, the coefficients k and r] may be ex- 
pressed as functions of k; they are evaluated for a specific value of k in 
Appendix C. Then, the last term of equation (16) is 

(~J(BA)\ ~* fTC^A _1 r~lBA)\ _1 ABA), ABA), ~JBA) 
\ jk ) \ J \ li J ks pl mn 

= [^4 BA) A BA) dS A) + 6 (tjJ { n BA) + ^ BA) ) + & nm df A) ] , (30) 

where £i = (rjiUx + r? 2 ^i + 2^2) , 6 = ?7i^2, and £ 3 = (^2 + ^1^2), with 
the coefficients 77 and uj defined in (29) and (18), respectively. 

5 Effective Moduli 

We use (6), (7), and (10) to write the increment in contact force as 



p (BA) = K iBA)^BA)^A) ^SA) + ^ 

where is given by (8) and (9) and 



k{BA) _ _ 9 ABA)-1 j(BA)rp 



(BA) 
j 



= -2 [r l [uJf A) d^d^ + co 2 (6 3m d[^ + 5 jn dS A) ) + uJ 2 S nm i 

- r 2 [*Jf A) d^d^ + K 2 (5 im $?*> + 5 jn d£ A) ) + ^nJf A) ] 

[ildf A) d^d[ B ^ + 6 {S jm ^ + 5 3 J r B n A) ) + ^nJf A) ] } E mn . 

With the increment contact force known, the increment in stress, 
may be calculated, making use of the identities 

%BA)*BA).~ 47T 



d^'d^dQ = f8 tJ 



and 
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Therefore, the incremental response of an isotropic, random aggregate of 
identical frictionless spheres is given by 



T 

= I 2 ! 1 - 2?TVi + 2^ 2 ) + 2^\k 1 + 2k 2 ) - 2^ 3 (6 + 2£ 2 )]^ 

+ [1 - 2ip~ 1 (uj 1 + 2uj 2 ) - l^- l u 2 + 2^~ 2 (/ti + 2k 2 ) + 10^~ 2 /t 3 
-2^" 3 (6 + 26) - 10^" 3 6] • 
/,From this, the effective moduli are 

ii E = ^K N {1 - 2 + 2cu 2 ) - ^- 2 ( Kl + 2/ca) + V" 3 (6 + 26)] } 

(31) 

and 

A B = — ^[l-2^- 1 (cu 1 + 7cu 2 ) + 2^- 2 (/ tl + 2/ t2 + 5/ t3 )-2V- 3 (6+26 + 56)], 
ona 

(32) 

where, we recall that ip = k/3 and the coefficients n, u, and £ are defined 
in (28), (18), and (30), respectively. The dependence of the effective moduli 
on the pressure occurs only through Kn and k, the other coefficients are 
functions only of the geometry. A dependence on p 1//3 enters through Kjy. 



6 Comparisons 

We compare the results of computer simulations with the predictions of the 
model. The numerical simulations were carried out using the code TRUBAL 
developed by Cundall (1988). Each simulation employed spheres of two dif- 
ferent radii: R\ = 0.105x 10~ 3 m and R 2 = 0.095 x 10~ 3 m in equal numbers. 
The shear modulus of the material of the spheres was /i = 2.9 x 10 10 Pa and 
the Poisson ratio was v = 0.2. This simulation employed a system of 10, 000 
spheres. The initial state was obtained in the manner described by Makse, et 
al. (1999). An initial random aggregate of frictionless spheres without con- 
tacts was homogeneously and isotropically contracted, bringing the spheres 
into contact, until a pressure p of 100 x 10 3 Pa was reached. In this state, 
k = 6.067 and v = 0.637. 

The numerical results for the shear and bulk moduli were 

fi E = 8 MPa and k e = 200 MPa. 

The predictions are based on identical spheres made of the same material 
with the average diameter d = 0.1995 x 10~ 3 m in an initial state with the 
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same coordination number, volume fraction, and confining pressure. Then, 
from equations (8) and (9), we first calculate the average normal stiffness: 



K N = 1.08 x 10 5 Pa-m. 

In this case, the effective moduli calculated from the average strain assump- 
tion are 

fj, E = \ E = 134 MPa, 

and 

n E = \ E + V = V = 223 MPa. 
3 P 3 P 

On the other hand, the effective moduli determined from equations (31) and 
(32) are 

fi E = -—K N (1 - 0.69) =41 MPa 
bird 

and 

\ E = -K N (1 + 0.57) = 208 MPa; 

bird 

so, the bulk modulus is 

k e = 235 MPa. 

The predicted value of the effective shear modulus is about 70% less than 
that resulting from the average strain assumption. This is encouraging. The 
increase in the bulk modulus over that based on the average strain assump- 
tion can be reversed by incorporating the anisotropic part of (A( BA )) 1 into 
the calculation. 



7 Conclusion 

We have considered a random aggregate of identical, frictionless, elastic 
spheres subjected to an initial confining pressure followed by a general incre- 
ment in strain. We have incorporated the deviation from the average of the 
difference in the displacement of a typical pair of particles into the contact 
force. We then determined the approximate value of this fluctuation in terms 
of statistical measures of the geometry of the packing and the contact stiff- 
nesses using force equilibrium in a way suggested by Jenkins (1997). We then 
made simple but plausible assumptions regarding the form of the simple and 
conditional probabilities that determined the geometry of the packing and 
used these to carry out the averages. Of particular interest was the quan- 
tity (A( Bj4 )) 1 J( BA ). This involved the product of averages and the average 
of the product of fluctuations. We found that the product of the averages 
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provides 46% of the reduction of the shear modulus and the average of the 
product provides the remaining 24%. Other sources of a reduction in moduli 
were not considered. 

For example, Trendadue (2001), apparently following Jenkins (1997), car- 
ries out a similar calculation of the bulk and shear moduli for elastic, fric- 
tional spheres and includes a reduction in moduli due to fluctuations during 
the initial isotropic compression and a reduction due to the retention of the 

anisotropic part of (A( Bj4 )) \ On the other hand, he does not include a con- 
tribution from the product of the averages. Our experience with the computer 
simulations indicates that the reduction of the shear modulus associated with 
the former two contributions are relatively small, but that associated with 
the the latter is significant. 

Independently, Paine (1997) introduced fluctuations in displacement and 
used the equilibrium equations to determine them in terms of the packing and 
stiffness. She then employed computer simulations, measured the statistical 
distributions, and calculated the reduction of the bulk modulus. 

We believe that the correct modeling of the statistical distributions is 
crucial to the prediction of the mechanical behavior of granular materials. In 
our modeling, we have assumed them to be the most homogeneous possible. 
This assumption is plausible and seems to be effective. We have focused on 
frictionless particles in an isotropic compressed state, because, for friction- 
less particles, the difference in the values of the shear modulus based on the 
average stain assumption and those measured in the numerical simulation 
is very large. Consequently, this material provides a good test of any im- 
proved theory. We believe that we have incorporated the essential features 
of the correct distributions into such a theory and have established an appro- 
priate basis for the extension of the theory to frictional particles and more 
complicated states of stress. 

Finally, it is natural to question why the average strain assumption pro- 
vides a relative good approximation to the bulk modulus and a relatively 
poor approximation to the shear modulus. An indication of why this is so 
can be obtained by considering a behavior of a simple model neighborhood of 
the sphere A in which sphere B is on the polar axis in the upper hemisphere 
and three other spheres in the lower hemisphere are distributed symmetri- 
cally about the polar axis at an angle = cos -1 (1/3) below the equator. 
Then, when this neighborhood is subjected to a isotropic compression, the 
relative displacement of the particles A and B is given by the average strain 
and there is no fluctuation. However, when the neighborhood is subjected 
to a deviatoric strain, there must be a deviation from the displacements as- 
sociated with the average strain, whose direction and magnitude varies with 
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the orientation of the neighborhood with respect to the principal axes of the 
strain, in order for the pair to be in equilibrium. That is, a typical neighbor- 
hood is closer to being equilibrated by a average isotropic strain than by an 
average deviatoric strain. 
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A Single Particle Relaxation 

Consider a specific particle, labelled A, which we take to be centered at 
the origin. It has contacts with particles centered at S- nA \ n = 1,2, ..N^ A \ 
Assuming that one of the particle centers is displaced by an increment u( nA ) 
the form of equation (6) is 

F> = K N [d k u K k 'Jdi , (33) 

where is given by 

K N = - d 1 

1 — v 

with 5 given by equation (9). 

As written, the total force on the specific particle, due to the sum of all 
the incremental contact forces is not zero: 

fiW ^ £ pip- A) Q 
n=l 

Accordingly, that particle will move to a new incremental position, X^. 
The generalization of equation (33) that takes into account the new position 
and orientation is 



p{nA) = K 



N 



1 {nA) 
(li , 



Now, the requirement that the particle is in equilibrium with its contact forces 
gives three linear equations in the three unknowns X^ A ^. It is straightforward 
to solve these equations numerically. 
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Having determined the new equilibrium position and orientation, one can 
show that the total increment in work done by the increment in contact forces 
on particle A is simply 



W {A) 



1 

2 



N (A) 



5 ( nA ) . (nA) 



U 



n=l 



(34) 



where X.^ is determined as described above. In order to calculate 

we make the average strain assumption that the displacement at the contact 

point is simply related to the macroscopic strain 



(nA) 



= E ij d J 



(nA) 



(35) 



Because we know the exact positions of each center, dS- nA \ from the simula- 
tions, we are able to evaluate equation (34) for each particle in the ensemble. 

We now evaluate the elastic moduli by setting the total deformation in 
the contacts equal to the macroscopic strain energy: 



1 N 1 
V ^ 2 



A=\ 



ii - -jl ) {E u ) 2 + 1ixE i3 E l3 



(36) 



where the sum is taken over all particles in the computational unit cell of 
volume V. The left-hand side of equation (36) can be evaluated once for a 
pure compression and once for a simple shear in order to deduce the values 
of R and p. The point of the exercise is to investigate the extent to which 
relaxation, at the single particle level, can explain the large reduction of the 
shear modulus relative to the prediction of the average strain approximation. 

If, in equations (34) and (35), we assume there is no relaxation (X^= 0), 
and if we replace the sum over contacts by an integral over a presumed 
uniform distribution of contact directions, we reproduce the average strain 
theory, equations (1) and (2), as detailed in Norris and Johnson (1997). 



B z a \p(at,P) 

Here, we determine the appropriate form of the distribution function z a \p(a, (3) 
for the intervals of polar angles 6 a and 6@ indicated in (25), given that particle 
B is on the pole. We assume that the distribution z a \p(a, j3) is independent 
of the circumferential angle <p a and that there are the same number of a con- 
tacts in each infinitesimal circumferential strip about the pole whether or not 
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the particle (3 intrudes on the strip. This is the simplest kind of homogeneity 
that applies in this case. 

We first consider the range of polar angles in the first integral: 

tt/3 <0p< 2tt/3 and vr/3 <0 a <0p + tt/3. 

In the subinterval ?r/3 < Op < vr/2, g((3) = [(k/2) — 1] /n — (k — 2)/2vr. For 
7r/3 < a < n I '2, the expected average number of a particles in the upper 
hemisphere, given the presence of (3 there, is (k/2) — 2, so 

/ z a \p(a,(3) sin a d<f> a d0 a 

J$>{e a ) 



7T 

= 2 ^ ^(0,) [tt - $(0 Q )] sin0 Q d0 a - 



If ^^(^a) [7r — 3>(0 Q )] is to be constant, then 

k — 4 

When 7r/2 < a < Op + ir/3, the average number k\ of a particles in this 
strip, given the presence of (3 in the upper hemisphere, is obtained from the 
proportion 

2tt r9/}+% 



k/2 : 2vr = k x : / / sin0 Q d0, 
Jo if 

as 

k 1 = - k - C o* (Vf). 

Then, 



--cos[0p + -)= / ^(a,/?) sin O a d(f) a dO a ; 

1 V 77 -/<I>fO 



and 



^a|/3(tt,/3) 



*(<? a ) 

k 



4(tt- $)' 

In the subinterval tt/2 < 0^ < 2tt/3, = (k/2)/2n = k/4n. For tt/3 < 
0« < 7r/2, the average number of a particles in the upper hemisphere, given 
the presence (3 in the lower hemisphere, is (k/2) — 1, so 

/ z al p(a,(3)smO a dO a d<p a: 
20 



and 



fc-2 
2(tt - ¥) 



For n / 2 < 9 a < 2n/3, the average number of a particles in this strip is 
(fc/4) - 1, so 



fc-4 



2ir-<S>(9 a ) 



^(cc,/?) sin 6 a d(j) a d8 a , 



and 



^a|/3(a,/3) 



fc-4 



Finally, for 2tt/3 < Q < ^ + tt/3, 



2tt /.Sfl 



SO 



fc/2 : 2tt = fc 2 



k 2 = -- 



' Jo h 



sin 9 a d9 a d(f> a , 



7Y 



COS [08 + - + 



fc 

2 



7T 



cos^ + - j + - 



/•27r-*(6> a ) 



and 



*(0a) 

A: 



^(a,/?) sm9 a d(f) a d0 a , 



4(tt-$)' 

The range of polar angles in second integral is 

vr/3 < ^ < 2tt/3 and ^ + vr/3 < Q < n. 

In the subinterval 7r/3 <9p< n/2, g{(5) = (k — 2)/27r; and in the subinterval 
tt/2 < ^ < 2tt/3, = fc/4vr. For ^ + tt/3 < 9 a < vr, given f3 in the first 
subinterval, the average number fc 3 of a particles in the strip is determined 
by 

2-7T P7V 



: 27r = fc 3 : / / sin 9 a d9 a d(f) a 
Jo Je^+f 



as 



Then 



fc 3 = 



1 + COS ( 6/3 + 



f^+o) =/ / 2a|/3(a,/3)sin0 a d0 Q d0 a , 



1 + COS 
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and 



k 

An' 



Then for 0p + 7r/3 < a < n and j3 in the second subinterval, the distribution 
function is the same. 

The range of polar angles in third integral is 

2tt/3 <0p<n and tt/3 < a < Op - tt/3. 

In the subinterval 27r/3 < Op < 5tt/6, g{(3) = k/An and the average number 
A; 4 of a particles in the strip 7r/3 < a < Op — tt/3 is determined by 



(k - 2) /2 : 7T = k 4 : 



2- 



sin O a dO a d(f) a 



as 



Then 



k 4 = (k- 2) 



- COS 00 



2tt 



z a |/?(a,/3) sin O a d(f) a dO a , 



and 



^a|/3(a,/5) = 



(fc-2) 
2tt 



In the subinterval 5vr/6 <O p <tc, 5r(/3) = k/4n. For vr/3 < Q < vr/2, the 
average number of a particles in the upper hemisphere, given the presence 
of (3 in the lower hemisphere, is (k/2) — 1, so 



z a] p(a,P) 



k-2 
2tt 



For 7r/2 < a < Op — 7r/3, the average number A; 5 of a particles in the strip 
is determined by 



k/2 :2n = k 5 : 



2w rOa-l 



sin O a dO a d(j) a 



as 



Then 



COS 



/ 7T\ f^"^ /" 27r 

yP~~%) = J J MP\ a iP) sinO a d(p a dO a , 
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and 

z a \p(a,/3) = — . 
The range of polar angles in the fourth integral is 

2tt/3 < Op < 7T and ^ - vr/3 < # a < 5vr/3 - 0^, 

where = k/4n. In the subinterval 0g — n/3 < 9 a < ir/2, with 27r/3 < 
Op < 5tt/6, 

(k - 2) /2 : tt = k 6 : / sin 6 a d9 a d(f) a , 



to Je fi -l 



k 6 = (k- 2) cos (00 - |) ; 
(fc - 2) cos I Op - -) = / / z^a, (3) sm 9 a d(p a dO a 



and ^ 

^a|/3(«,/5) 



2(tt-$)' 
Similarly, from 



r>2vr 

k/2 : 2tt = k 7 : I I sin 9 a d9 a d(f) 



so 



and 



^ = —008 ^-61^; 
k 7 = / z a |0(a,/3)sin0 a d0 a d0 a 



z, 



4(tt - $) 

Finally, for 9p — n/3 < 9 a < 5^/3 — 0^, with 5n/6 < 9p < it, the distribution 
function is the same. 



23 



C Details of the Calculation 



C.l A^)'J(^)' 

We consider first the integrals in equation (25) in which the limits of the 
integration over <p a depend upon $ and define 

Qfmn = / (3j(3id<f)(3 / aia m a n d(t) a . 

Jo J$+<t>f] 

When expressed in terms of the angles, this is 



q{BA) = 
jmn 

7T [2(tt - $) sin 2 Q cos 9 a cos 2 dp - sin $ sin 26 p sin 3 6 a ] 5 mn df A) 
+ - [2(n - $) sin 2 a cos a sin 2 0^ - sin 2$ sin 2 9 a cos 9 a sin 2 0^ 

-2 sin $ sin a cos 2 6 a sin 20^] ^ mj rff A) + 5 jn d { * A) ^ 
+ - [4 sin $ sin 0^ cos dp sin 3 a + 8(7r - $) cos 3 a cos 2 0/3 
-4(tt - $) sin 2 a cos a + 2 sin 2$ sin 2 a cos a sin 2 0^] d[^df A) d^ A) . 
With this, the integral 

R{ j^n = I I g(P)z a \p( y a,p)p j p i a i a m a n dn a dnp 



in equation (25) can be written as 



3 



0/3 + 1 



= / / g(P)z a \p(a, P)Q ( S£ sin a sin 6pd9pdO a 
+ / / / / g{l3)z a \p(a ) f3)f3 j f3 i a i a m a n dVt a dVtp 

J\ JOg + ^JO JO 



73+3 

7T 

3 



^27r /*27r 

+ 11 I g{f3)z a \p(a ) f3)f3 j f3ia i a m a n dVL a dVtp 



'27T J 7T 

3 3 

5 77 



'0 ./0 



+ / c/(/?)z a |^(a,/3)Q^ ) sin0 a sin0 /3 rf0 / 3ci0 a , 

We use the appropriate distribution function for each range of variation of 
the polar angles and distinguish the integrals that involve Q( BA ) from those 
not. 
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Then, for example, 



t\p(a, P)Qfr£! sin 9 a sin 9gd9gd9 a 



r 2 -i rh+l 
/ / 9(P)z, 

= fc^(jb_4) [ 2 [ 2 Qjmn ^ sing, sin 9 a d9gd9 t 

3 3 

- — Jmi 1 sin 9g sin 9 a d9gd9 a 
i [it — $) 




+ 



-(.-2)^ ^ ^sm^sm W ^ 

_2tt „2tt ^(BA 



'a 



/ - — Jm " . sin flg sin 9 a d9gd9 a 



^ i67r /| y^r 



/ - — j tyiti sm9gsm9 a d9gd9 a . 
1 2* (ir-$>) 1 p 



2 3 v ' 

evaluation the integrals over 9g and # a numerically, we obtain 
I I g(P)z al p(a,/3)Q^ sin9 a sm9gd9gd9 a 

' in^(k- 4) + 0.57^> -0.38^^ 
4tt V ' 8tt 8tt 



+0.13— (A; -4) + 0.19— 
16tt ; 167T 



+ 



Q n (fc_2) _ 4) _ Q 23 (fc_2) k + Q13 k _ 2 ) 

47T 87T 57T 

jU n -r u jn u m I 



-0.11 A (t _ 4) -0.11^ 



+ 



Q n (fc_2) _ 4) _ Q lg (fc_2) fc + ^ _ 2) 

47T 87T 87T 



+ 0.ll-^-(jfe-4) + 0.01— 

16tt y 167T 



5 d {BA) 

u mn u 'j 
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In the same way, 



if- y /3 



I 2n . .. 




g((3)z a \p(a, P)Q^ sin 9 a sin 9pd9 d9 a 



k(k - 2) tl Q) 



57T /2tt /a _tt_ (7T 
3 



,(BA) 



" ]mn ^ sin 9 a sin 9pd9pd9 a 



k 2 fir ff-^ q(^) 
+ — — / / 7 sin 9 a sin 9pd9pd9 a 

l07T J^r J| (7T — 4>) 



k 2 r r¥- e Q {B J^ 

+ — — / / / sin 9 a sin 9 p d9pd9 a 

16tt J ¥ ^ ? (tt - $) 



fc(fc-2) 



8tt 

+0.105 mn dT A) 



+ 



k 



(-0.234^)^^) - 0.055 mn rff A) ) 



167T 

In the integrals that are independent of $, we can easily carry out the 
integration of 

P jm^ = y PjPid^PT^. J aia m a n d(j) a : 

^ = d£ A) d*f A) di BA) [cos 2 ^ (cos 3 Q - ? cos 9 a sin 2 6> a ) 
— - sin 2 0p cos # Q sin 2 9 a + cos 2 0g cos 9 a sin 2 a 



- cos 2 9/3 cos Q sin 2 9 a df A) 5 mn 



l - sin 2 cos 9 a sin 2 Q (^d?^ + 6 nj d£ A ^ 
With this, 



+ 4 



. 2tt 
3 



2tt ^2tt 

g((3)z a \p(a, (3)(3j(3iaia m a n dQ a dQp 




Jep+i Jo Jo 
M^)/ f r P fJ^9 a sin9pd9pd9p 



3 " W P^3 



+ T / 3 P i-n ) sin a sin 

4 J| Je 0+ i 



26 



and 



k(k-2) 
' 8^ 
-0.03<& BA) 6. 



OlOd^df^d^ - 0.09 (5 jm d^ + 5 jn d£ A) ) 



+ 



16tt 



omd^d^d^ - o.oi (Vf 4 ) + sJ? A) ) 



3 P 



7r p2tv p2tt 



./O 



+ 



k{k-2 



2lT I 7T 

3 3 



2 / 52L / 

"63 



2 /*7T 



4 / 57T / 7T 

'62 



+o.o8dr ; ^ m „ 



+ 



fc(ife-2) 
8^ 

k{k-2) 
8^ 



-0.12^)^4^ + 0.03 fe m ^ + * iB <g"> 



-0.07d^df A) dl BA ^ + 0.01 + 5 jn d£ A > 



+0.08dT A> 8 mn 



+ 



16tt 



Omd^df^d^ - 0.02df A) 6 mn 



Then 



167Ti? 



(BA) 
jmn 



-- - [0.52(fc - 2)(ife - 4) + 0.10fc(fc - 2) 
-0.13jfe(ifc - 4) - 0.01A; 2 ] df A) d^d^ 
+ [0.44(A: - 2)(ife - 4) - 0.24A;(A; - 2) 

-0.11A;(A; - 4) - 0.14A: 2 ] (6 jm S^ + M£f A) ) 
- [0.44(fc - 2)(ib - 4) - 0.42A:(fc - 2) 
-0.11Jfe(A; - 4) + 0.04A; 2 ] 5 mn df A) 
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Finally, we write the first integral in (20) as 



Sjmt^= / / F(a,(3)(3 j (3ia i a rn a n dQ Q dQf 3 



, R (BA) 



C.2 A(^)'A^)' 



For the calculation of A^^'A^f A ^ we introduce the integral 



= / (3j(3 s d<j)p / a s aid(j) a . 

JO J$+(I>I3 



When expressed in terms of the angles, 

Fji — —- [sin 2$ sin 2 a sin 2 Op + sin $ sin 20 a sin 20 p 
-2(7r-$)sin 2 a sin 2 Op] S jt 
+ - [sin 2$ sin 2 9 a sin 2 Op - sin $ sin 29 a sin 29 p 
-2(tt - $) sin 2 6> Q sin 2 6^ + 8(vr - $) cos 2 a cos 2 6*3] dj-d, 
With this, the integral 

Y ji BA) = / F!jf A) g(f3)z a \p(a,{3){3 j f3 s a s a m dtt a dttp 
in equation (25) can be written as 

Y^ A) = J n ^ F^ A) g({3)z al p(a,P)sm0 a sm0pd0pd0 a 



3 3 
2tt 
3 



' 1 2-7T /*2"7T 

+ 11 I g(/3)z a \p(a, p)pjp s a s aidn a dflp 



3 "^ + f^0 Jo 
7T ^e fl -? /.2tt /.2vr 



+ y 2 J w J J g{l3)z a \p{a ) f3)f3 j f3 s a s aidVL a dVtp 
+ / / F^f A) g((3)z a \p(a,(3)sin0 a sm0pd0pd0 a , 
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Then, for example, the first integral is 



I / g((3)z a \p(a,f3)Ff A) sm6 a sm6pd6 p d6 a 



~^-i k - 4 ) / 2 / 2 7 ~ ^ Sil1 ^ Sin 



3 3 



(fc-2) 




+ — -A; / / — ^ sin ^ sin 9 a d9pd9 a 



87T 7| 7| (7T - $) 

+ —(k-2) / iZ — sm6f3sm6 a d6f3d6 a 

o7T 7| 7| (7T — <3? J 

lb /•¥/•¥ F ( f A) 

+ —(k-4) / J ' sm9p sin 9 a d9pd0 a 

l07T (7T — <P) 

+ — / / k\ sm9 p sm9 a d9pd9 a . 

107T ,/| J^r (7T — <P) 

Upon evaluating the integrals over # Q and 0p numerically, we obtain 
/ /5)i^7 sin 6 a sin OpdOpdOa 

A/3 

{k ~ 2 \k-A) (0.495,, - 0.544 BA) ^ A) ) 



47T 

+ 



^— — U fo.71*,-, - O.SQd^d^ 
87T V J J 

+ ^ (Jb _ 2 ) (0.555,, - 0.47df A) rff A) ) 
+ (0.49^-0.544^^)) 
+ ^(0.17^-0.16^^). 
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Similarly, for the second integral, 



( Ft A \ 

9{P)z a \p{a, (3)F\ sin 6p sin 6 a d9 a d9p 

! 2ir In tt J 
— J( >0-3 

k(k-2) H Ff A hin6psmO a d6 a d6p 




8tt J 2, J^_, (tt-$) 



roiT/b p- 



167T J27V, Jn {jl — $) 



+ 



3 2 



6 P 3 



(0.16^-0.124 BA) rff A) ) 



fc(ife-2) 
8^ 

16tt "I ) 



+ — fo.14^, + o.o8$ SA) d, (iM)> ) 



0.0% -0.14^^). 



In order to calculate the last two integrals, we first introduce 



1*271 /»27r 

Cjf A) = — J pjpsdcpp— J a s aid(j) a 



- sin 2 6/3 sin 2 9 a 8ji + (cos 2 6*^ cos 2 6 a — — sin 2 9p sin 2 # Q ] d^ A ^d[ BA 



Then, 



/ / / / g(f3)z al/3 (a,p)p j p s a s aidQ a dQp 
J% Jdp+j Jo Jo 



— — — ^ / C\f A) Sin6 a sin6pd6 p d6p 



2 

3 

2tt 




+ T" / 3 / C< 7 ( f A) Sil1 ^ Sil1 

= ^ (0.1^ + 0.01^>c?-)) + ^0.01^, 
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and 



g((3)z a \ P (a, (3)Cf A) sin 8 P sin 8 a d8 a d8p 



2_7T / 7T 

3 3 



k(k-2) 



2-7T / 7T 

3 3 



C-f A ' ) sin 8/3 sin 6 a d6 a d6 a 



+ 



+ 



k(k-2) r 



2 

,2 



57T / 7T 

6 3 



Cjf A ^ sin 6?g sin 8 a d8 a d8p 



k 2 r r 

4 J 7T 



Cjf ^ sin 0g sin 8 a d8 a d8p 



87T V J / 

+ Miz^l ( .02^ + 0.04^^>) + ^0.01*,-,. 

So, 

= [1.9Q(k - 2){k - 4) + 3.30A;(A; - 2) + 0.49A;(A; - 4) + 0.32A; 2 ]^ 7 

- [2.16(ib - 2){k - 4) + 2.30A;(A; - 2) + 0.54A;(A; - 4) - 0.06A; 2 ] df A) d\ BA) . 

(39) 

Finally, we write the first integral in (23) as 

H jf A) = [ [ FfaP^jPiOiaidnccKlfi (40) 

= r/r ) + «i^ + 5 2 C ) rfr ) . 

C.3 Final evaluation 

In order to compare the values of the elastic moduli obtained by numerical 
simulation with what predicted by the theory, we have to evaluate all the 
previous quantities. These are functions of the coordination number, which 
is assumed to be k = 6.07. Therefore, with equation (37), we can write 

RfrrS = X,d { i A) df A) d^ + X2 (5 ml $ n BA) + M^) + X^ndf A \ 
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where Xi — —0.10, X2 = —0.17, and % 3 = 0.13 and, with equation (39), 

Y^= Pl 5 jl + P2 df A ^ BA \ 

where p\ = 2. 31 and p 2 = —1- 58. 

We can also evaluate all of the coefficients involved in the tensor formulas 
introduced earlier. So we have a\ = 1.94, a 2 = 0.24, lo\ = 0.78, uj 2 = —0.16, 
if> = 2. 02, a 2 = -0.76, and Zj x = -0.22. 

Using equations (38), (37), and (21), we derive 

SfJ = aj^d^d^ + a 2 (5 m3 dl BA) + Sj { ° A) ) + a,5 m Jf A ^ (41) 

where a x = —0.32, a 2 = —0.33, and a 3 = —0.0 3 and 

H^^b^ + bjf^d^, (42) 
where 61 = 4. 25 and b 2 = —2. 34. Then, from equation (26), 



ABA), j{BA), 



ji imn 



K 2 



where K\ = —0.31, k 2 = —0.02, and K3 = 0.15. This permits the calculation 
of the second term in equation (16). In a similar way, from equation (27), 



) 

where 7/1 = 0.49 and r\ 2 = 0.0 3. Then, with this and equation (18), the last 
term in equation (16) is 



(BA), A (BA)l T (BA) 

(BA) 



(^r)"(^r) _1 (4f 4, )"4f"4f 4 "^ 

where £1 = 0.40, £2 = —0.08, and £3 = —0.0 8. Finally, the first term in 
equation (16) is the simple average 



(^Y^n = r l [^df A ^d^ 

+u 2 (6 jn dg A > + 8 3 J^ + 5 n Jf A) ) 

that follows from equations (19) and (18). 
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